#include<iostream>
#include<sys/time.h>
#include<arm_neon.h>
using namespace std;
const int N=5001;
float m[N][N];
void m_reset(int n){
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++)
            m[i][j]=0;
        m[i][i]=1.0;
        for(int j=i+1;j<n;j++)
            m[i][j]=rand();
    }
    for(int k=0;k<n;k++)
        for(int i=k+1;i<n;i++)
            for(int j=0;j<n;j++)
                m[i][j]+=m[k][j];
}

//Neon第一处并行
//优化第一个内嵌循环中的A[k, j] := A[k, j]/A[k, k]
void neon1(int n) {
    struct timeval h, t;
    float time = 0.0;
        gettimeofday(&h, NULL);
        for(int k = 0; k < n; k++) {
            float* va = new float[4];
            *va = *(va+1) = *(va+2) = *(va+3) = m[k][k];
            float32x4_t vt = vld1q_f32(va);
            //将四个单精度浮点数从内存加载到向量寄存器
            delete[] va;
            for(int j = k + 1; j + 3 < n; j += 4) {
                float32x4_t va = vld1q_f32(&m[k][j]);
                va = vdivq_f32(va, vt);//向量化除法运算
                vst1q_f32(&m[k][j], va);//储存到内存
            }
            for(int j = n - n%4; j < n; j++) {
                m[k][j] = m[k][j]/m[k][k];
            }//结尾几个元素串行除法运算
            m[k][k] = 1.0;
            for(int i = k + 1; i < n; i++) {
                for(int j = k + 1; j < n; j++) {
                    m[i][j] = m[i][j] - m[i][k]*m[k][j];
                }
                m[i][k] = 0;
            }
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    cout<<"neon1: "<<time<<" ms"<<endl;
}
//Neon第二处并行
//优化双层 for 循环中的A[i, j] := A[i, j] − A[i, k]×A[k, j]
void neon2(int n) {
    struct timeval h, t;
    float time = 0.0;
        gettimeofday(&h, NULL);
        for(int k = 0; k < n; k++) {
            for(int j = k + 1; j < n; j++) {
                m[k][j] = m[k][j]/m[k][k];
            }
            m[k][k] = 1.0;
            for(int i = k + 1; i < n; i++) {
                float* aik = new float[4];
                *aik = *(aik+1) = *(aik+2) = *(aik+3) = m[i][k];
                float32x4_t vaik = vld1q_f32(aik);//将四个单精度浮点数从内存加载到向量寄存器
                delete[] aik;
                for(int j = k + 1; j + 3 < n; j += 4) {
                    float32x4_t vaij = vld1q_f32(&m[i][j]);
                    float32x4_t vakj = vld1q_f32(&m[k][j]);
                    float32x4_t vx = vmulq_f32(vaik, vakj);//向量化乘法运算
                    vaij = vsubq_f32(vaij, vx);//向量化减法运算
                    vst1q_f32(&m[i][j], vaij);//存储到内存
                }
                for(int j = n - n%4; j < n; j++) {
                    m[i][j] = m[i][j] - m[i][k]*m[k][j];
                }//结尾几个元素串行计算
                m[i][k] = 0;
            }
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    cout<<"neon2: "<<time<<" ms"<<endl;
}

//Neon全部并行（不对齐）
void neon_normal(int n) {
    struct timeval h, t;
    float time = 0.0;
        gettimeofday(&h, NULL);
        for(int k = 0; k < n; k++) {
            float* va = new float[4];
            *va = *(va+1) = *(va+2) = *(va+3) = m[k][k];
            float32x4_t vt = vld1q_f32(va);//将四个单精度浮点数从内存加载到向量寄存器
            delete[] va;
            for(int j = k + 1; j + 3 < n; j += 4) {
                float32x4_t va = vld1q_f32(&m[k][j]);
                va = vdivq_f32(va, vt);//向量化除法运算
                vst1q_f32(&m[k][j], va);//储存到内存
            }
            for(int j = n - n%4; j < n; j++) {
                m[k][j] = m[k][j]/m[k][k];
            }//结尾几个元素串行计算
            m[k][k] = 1.0;
            for(int i = k + 1; i < n; i++) {
                float* aik = new float[4];
                *aik = *(aik+1) = *(aik+2) = *(aik+3) = m[i][k];
                float32x4_t vaik = vld1q_f32(aik);//将四个单精度浮点数从内存加载到向量寄存器
                delete[] aik;
                for(int j = k + 1; j + 3 < n; j += 4) {
                    float32x4_t vaij = vld1q_f32(&m[i][j]);
                    float32x4_t vakj = vld1q_f32(&m[k][j]);
                    float32x4_t vx = vmulq_f32(vaik, vakj);//向量化乘法运算
                    vaij = vsubq_f32(vaij, vx);//向量化减法运算
                    vst1q_f32(&m[i][j], vaij);//存储到内存
                }
                for(int j = n - n%4; j < n; j++) {
                    m[i][j] = m[i][j] - m[i][k]*m[k][j];
                }//结尾几个元素串行计算
                m[i][k] = 0;
            }
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    cout<<"neon_normal: "<<time<<" ms"<<endl;
}

//Neon全部并行（对齐）
void neon_algn(int n) {
    struct timeval h, t;
    float time = 0.0;
        gettimeofday(&h, NULL);
        for(int k = 0; k < n; k++) {
            float* va = new float[4];
            *va = *(va+1) = *(va+2) = *(va+3) = m[k][k];
            float32x4_t vt = vld1q_f32(va);//将四个单精度浮点数从内存加载到向量寄存器
            delete[] va;
            for(int j = k + 1; j < k + 4 - k%4; j++) {
                m[k][j] = m[k][j]/m[k][k];
            }//串行计算至对齐
            for(int j = k + 4 - k%4; j + 3 < n; j += 4) {
                float32x4_t va = vld1q_f32(&m[k][j]);
                va = vdivq_f32(va, vt);//向量化除法运算
                vst1q_f32(&m[k][j], va);//储存到内存
            }
            for(int j = n - n%4; j < n; j++) {
                m[k][j] = m[k][j]/m[k][k];
            }//结尾几个元素串行计算
            m[k][k] = 1.0;
            for(int i = k + 1; i < n; i++) {
                float* aik = new float[4];
                *aik = *(aik+1) = *(aik+2) = *(aik+3) = m[i][k];
                float32x4_t vaik = vld1q_f32(aik);//将四个单精度浮点数从内存加载到向量寄存器
                delete[] aik;
                for(int j = k + 1; j < k + 4 - k%4; j++) {
                    m[i][j] = m[i][j] - m[i][k]*m[k][j];
                }//串行计算至对齐
                for(int j = k + 4 - k%4; j + 3 < n; j += 4) {
                    float32x4_t vaij = vld1q_f32(&m[i][j]);
                    float32x4_t vakj = vld1q_f32(&m[k][j]);
                    float32x4_t vx = vmulq_f32(vaik, vakj);//向量化乘法运算
                    vaij = vsubq_f32(vaij, vx);//向量化减法运算
                    vst1q_f32(&m[i][j], vaij);//存储到内存
                }
                for(int j = n - n%4; j < n; j++) {
                    m[i][j] = m[i][j] - m[i][k]*m[k][j];
                }//结尾几个元素串行计算
                m[i][k] = 0;
            }
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    cout<<"neon_algn: "<<time<<" ms"<<endl;
}

int main(){
    int n=64;
    m_reset(n);
    neon1(n);
    m_reset(n);
    neon2(n);
    m_reset(n);
    neon_normal(n);
    m_reset(n);
    neon_algn(n);
    return 0;
}